%% SOLVES ODE SYSTEM IN LIMITED-COMMITMENT CONTRACT
function [V, Z, C, Beta, K, Eta, f, ExRet, mu_v, vol_phi, mu_k, mu_c, vol_k, vol_c, Suff] = limited_system_solver(params, grid_params, num_params)
%% PARAMTERS
temp = num2cell(params);
[rho, r, ra, alpha, sigma, lambda, M, eta_min, eta_max] = temp{:};

[n_grid, ~, PHI_mat, ~, dphi] = grid_params{:};

[max_iter, D, thr] = num_params{:};

PHI = PHI_mat(1,:)';

%% Initialize matrices

V_out = (((ra-r)*(1-rho)/rho)^(rho/(1-rho))); %Value for shutting down the firm

V00 = V_out*ones(n_grid, 1);

C0 = 0.01*ones(n_grid,1);
Beta0 = 0.01*ones(n_grid, 1);
Eta0 = eta_max*ones(n_grid, 1);

%% derive implicit information rent

Z0 = limited_inforent_solver(C0, Beta0, Eta0, params, grid_params);


%% ITERATIONS

iter = 1;
reiterate = 1;
while reiterate
            
    tic        

    %% Find optimal policy.    
    [C, Beta, Eta, K, f, ExRet, mu_v, vol_phi] = limited_optimal_policy(V00, Z0, C0, Beta0, Eta0, params, grid_params);
    
    %% Find info rent
    Z = limited_inforent_solver(C, Beta, Eta, params, grid_params);
    
    %% Find value function
    V = limited_hjb_solver(V00, Beta, f, mu_v, vol_phi, params, grid_params, num_params);
   
    %% Check change
    vmin = min(V);
    vmax = max(V);
    
    [posmin] = find(V == vmin,1); 
    [posmax] = find(V == vmax,1); 
    
    disp(['Min of ', num2str(vmin), ' at ', num2str(posmin)])
    disp(['Max of ', num2str(vmax), ' at ', num2str(posmax)])
    
    disp(['Step size of ', num2str(D)])

    changemat = abs(V./V00 -1);
    test = max(changemat);
    [poschange] = find(changemat == test,1);
    disp(['Change of ', num2str(test,'%.4e'), ' at ', num2str(poschange)])
    
    
    if  test < thr
        reiterate = 0;
        disp(['Converged at iteration ',num2str(iter)])
    elseif sum(isnan(V))>0
        disp('NaN appeared');
        break
    else
        V00 = V;
        disp(['Iteration ',num2str(iter)])
        C0 = C;
        Eta0 = Eta;
        Beta0 = Beta;
        Z0 = Z;
    end
       
    iter = iter+1;
    
    toc
    
    disp(' ')
    disp(' ')
    
end

%% Dynamics of Capital

Kp = K;

%Initialize matrices
Kphi_f = zeros(n_grid,1);
Kphi_b = zeros(n_grid,1);
Kphi_c = zeros(n_grid,1);

% First Differences w.r.t. phi
Kphi_f(1:end-1) = (Kp(2:end) - Kp(1:end-1))/dphi;
Kphi_b(2:end) = (Kp(2:end) - Kp(1:end-1))/dphi;
Kphi_c(2:end-1) = 0.5*(Kphi_f(2:end-1) + Kphi_b(2:end-1));

%Second Differences w.r.t. phi
Kphiphi = (Kphi_f - Kphi_b)/dphi;

% DRIFT
mu_k = 0.5*Kphiphi.*vol_phi.^2;

% VOLATILITY
vol_k = Kphi_c.*vol_phi;

%% Dynamics of Consumption

Cp = C;

%Initialize matrices
Cphi_f = zeros(n_grid,1);
Cphi_b = zeros(n_grid,1);
Cphi_c = zeros(n_grid,1);

% First Differences w.r.t. phi
Cphi_f(1:end-1) = (Cp(2:end) - Cp(1:end-1))/dphi;
Cphi_b(2:end) = (Cp(2:end) - Cp(1:end-1))/dphi;
Cphi_c(2:end-1) = 0.5*(Cphi_f(2:end-1) + Cphi_b(2:end-1));

%Second Differences w.r.t. phi
Cphiphi = (Cphi_f - Cphi_b)/dphi;

% DRIFT
mu_c = 0.5*Cphiphi.*vol_phi.^2;

% VOLATILITY
vol_c = Cphi_c.*vol_phi;

%% SUFFICIENT CONDITION
% Holds if Suff matrix is positive

Zphi = 0.5*( Z - [0;Z(1:end-1)] + [Z(2:end);0] - Z )/dphi;
vol_Z = Eta.*PHI.*(1-PHI).*Zphi;

scaled_omega = (1-rho)*Beta.*Z + vol_Z;
Suff = scaled_omega - Eta.*(1-2.*PHI).*Z;



end
